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ABSTRACT 

The abundance of galaxy clusters can constrain both the geometry and growth of 
structure in our Universe. However, this probe could be significantly complicated by 
recent claims of nonuniversality - non-trivial dependences with respect to the cos¬ 
mological model and redshift. In this work we analyse the dependance of the mass 
function on the way haloes are identified and establish if this can cause departures 
from universality. In order to explore this dependance, we use a set of different N-body 
cosmological simulations (Le SBARBINE simulations), with the latest cosmological 
parameters from the Planck collaboration; this first suite of simulations is followed by 
a lower resolution set, carried out with different cosmological parameters. We identify 
dark matter haloes using a Spherical Overdensity algorithm with varying overdensity 
thresholds (virial, 2000pc5 lOOOpc, SOOpc^ 200pc and 200p5) at all redshifts. We notice 
that, when expressed in term of the rescaled variable z/, the mass function for virial 
haloes is a nearly universal as a function of redshift and cosmology, while this is clearly 
not the case for the other overdensities we considered. We provide fitting functions for 
the halo mass function parameters as a function of overdensity, that allow to predict, 
to within a few percent accuracy, the halo mass function for a wide range of halo 
definitions, redshifts and cosmological models. We then show how the departures from 
universality associated with other halo definitions can be derived by combining the 
universality of the virial definition with the expected shape of the density profile of 
halos. 
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1 INTRODUCTION 

In the Cold Dark Matter (CDM) model, structures - up 
to protogalactic scales - form through the amplification of 
small density fluctuations via gravitational instability (Lon- 
gair 1998; Springel et al. 2005; Mo et al. 2010; Angulo et al. 
2012). Dark matter haloes are objects which have been able 
to break away from the expanding background, and collapse 
(Press & Schechter 1974; Springel et al. 2001b). Small haloes 
form first, before merging with one another to form ever 
more massive ones in a hierarchical process (Lacey & Cole 
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1993, 1994). As a result of repeated mergers, dark matter 
haloes grow more massive in time (Tormen et al. 2004). 
Haloes hosting galaxy clusters represent the most massive 
and recently formed structures in our Universe. The forma¬ 
tion and the merger rates of dark matter haloes are sensitive 
to the expansion history of the Universe and so can be used 
to constrain cosmological parameters (Lacey & Cole 1993, 
1994; Moreno et al. 2008). 

In particular, different theoretical studies have shown 
that the halo mass function and its evolution are important 
probes of the very early Universe, its expansion history, and 
the nature of gravity (Press & Schechter 1974; Bond et al. 
1991; Lacey & Cole 1993; Sheth et al. 2001). These have 
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shown that, in appropriately scaled units, the mass function 
can be written in an approximately universal form which 
is independent of power spectrum and expansion history. 
Although this universality is only expected to be approxi¬ 
mate (Musso & Sheth 2012; Paranjape et ah 2013), it vastly 
simplifies the process of constraining cosmological param¬ 
eters from observational datasets, so it has served as the 
basis for fitting functions whose parameters are calibrated 
using numerical simulations of the dark matter (Sheth & 
Tormen 1999; Jenkins et al. 2001; Warren et al. 2006; Sheth 
& Tormen 2002; Del Popolo & Gambera 1998, 1999). As 
simulated datasets have grown, it has become possible to 
quantify (small) departures from universality (Tinker et al. 
2008; Crocce et al. 2010; Manera et al. 2010; Wu et al. 2010; 
Courtin et al. 2011; Corasaniti & Achitouv 2011; Murray 
et al. 2013; Watson et al. 2013). Several recent works have 
also concentrated on the small but significant modifications 
by neutrinos (Castorina et al. 2014), coupling between dark 
matter-dark energy (Cui et al. 2012; Giocoli et al. 2013) 
and different baryon physics (Gui et al. 2014; Bocquet et al. 
2015). 

Many cosmological constraints have been obtained from 
cluster counts (Vikhlinin et al. 2009; Rozo et al. 2010; Planck 
Gollaboration et al. 2013a) and other galaxy cluster prop¬ 
erties (Evrard et al. 2008; Ettori et al. 2009; Giocoli et al. 
2012c). Eor an extensive review, see Borgani & Kravtsov 
(2011). However in the near future, many wide-held surveys 
are expected to use the cluster mass function to constrain 
cosmological parameters (Pillepich et al. 2012; Laureijs et al. 
2011; Sartoris et al. 2015; Boldrin et al. 2015). In the light 
of these, a mass function calibrated to an accuracy of a few 
percent, and hexible enough to account for different halo 
identihcation dehnitions, is of primary importance. 

Departures of the halo mass function from universality 
may depend on how haloes are dehned. One of the main 
goals of this paper is to explore this dependance. Our ma¬ 
jor result is that, if haloes are dehned using the virial den¬ 
sity, and the htting formula includes a parameter which is 
related to this - as was done by Sheth & Tormen (1999) 
- then the mass function can be considered universal to 
within a few percent. In this respect, our hndings conhrm 
those of Gourtin et al. (2011): departures from universality 
result from ignoring the redshift and cosmology dependance 
of these quantities. Moreover, the departures from univer¬ 
sality associated with other halo dehnitions can be derived 
from combining the universality of the virial dehnition with 
knowledge of the enclosed density prohle around haloes. 

The paper is organised as follows: In Sec. 2 we describe 
the cosmological simulations we use for our study. It also de¬ 
scribes our halo hnder. We present our reference model for 
the halo mass function in Sec. 3. In Secs. 4 and 5 we discuss 
the universality associated with using the virial overdensity 
to dehne haloes and show how the parameters of the halo 
mass function depend on how haloes are dehned. The bulk 
of this analysis is for spherical halos: Appendix A shows how 
our results are modihed if halos are allowed to be ellipsoidal, 
and Appendix B describes how a number of more technical 
details affect our measurements. In Sec. 6 we show how the 
nonvirial halo mass functions can be derived from combin¬ 
ing knowledge of the dark matter density prohle with the 
virial mass function. In Sec. 7 we present some comparisons 
with previous works. We discuss our results and conclude in 


Sec. 8. Our analysis suggests that the most massive end of 
the halo mass function is particularly simple, in the sense 
that it can be described by a function with fewer free pa¬ 
rameters. Appendix G describes how this impacts cluster 
cosmology. All logarithms where not explicitly stated in the 
text are in base ten. 

2 THE NUMERICAL SIMULATIONS 
2.1 Le SBARBINE simulations 

Le SBARBINE simulations are a set of six dark-matter- 
only cosmological simulations run by the Padova cos¬ 
mology group. These simulations follow the evolution of 
1024^ particles, whose motions are assumed to be driven 
by gravitational instability, using the publicly available 
code GADGET-2 (Springel 2005). The assumed background 
cosmology and initial conditions for these runs are con¬ 
sistent with recent Planck results (Planck Gollaboration 
et al. 2013b) (hereafter Planckl3). In particular we have 
set: Drn = 0.307, Qa = 0.693, as — 0.829 and H = 
100/ikm s~^ h~^Mpc with h = 0.677. The initial power 
spectrum was generated using the GAMB code (Lewis 
et al. 2000), and initial conditions were produced by per¬ 
turbing a glass distribution with N-GenIG (http://www. 
mpa-garching.mpg.de/gadget); the realisations have been 
carefully chosen, in order to follow the initial power spec¬ 
trum even at large scales and thus to reduce the differences 
between the linear spectrum and that measured from the 
simulations. 

The parameters of our main simulation set are listed in 
Table 1. We used a different seed for the random number 
generator which sets the initial conditions of each simula¬ 
tion, so as to have a sample of independent realisations. Al¬ 
though each box contains the same number of dark matter 
particles (1024^), the comoving box lengths are different, so 
the mass resolution in each box is different. The box sizes 
were chosen so that the set provides good mass resolution 
down to 10^ h~^MQ. 

Le SBARBINE simulations are complemented by a set 
of lower resolution runs having different cosmological pa¬ 
rameters. These all have 512^ dark matter particles. In par¬ 
ticular, for each set of cosmological parameters we ran two 
simulations: one with box-size 150/i“^Mpc and another with 
1000/i“^Mpc. These were chosen to ensure good resolution 
both for intermediate and high mass haloes. This lower res¬ 
olution set was produced specifically to test the universality 
of the halo mass function with respect to the cosmological 
model (see Section 4.2). The parameters of these other sim¬ 
ulations are listed in Table 2. 

In addition, we also re-ran three simulations for which 
the initial conditions were generated using a second or¬ 
der Lagrangian Perturbation Theory (2LPT) algorithm ^ 
(Grocce et al. 2006). These are copies of Bice (1024^ parti¬ 
cles) and of the two 512^ simulations with the WMAP7 cos¬ 
mology (Komatsu et al. 2011) namely (wmap7 and wmap7- 
big). As we discuss in Appendix B4, there are small differ¬ 
ences - not exceeding 5% - between the mass functions in 
simulations with Zel’dovich versus 2LPT initial conditions. 

^ http://cosmo.nyu.edu/roman/2LPT 
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Main set of simulations 

name 

box [ h ^Mpc] 


mp[MQh 

soft [kpc h 

1 

o 

77 

II 

o 

Nh>30o{z — 0) 

colour 

Ada 

62.5 

124 

1.94 X 10^ 

1.5 

2264847 

103852 

green 

Bice 

125 

99 

1.55 X 10® 

3 

2750411 

129674 

cyan 

Cloe 

250 

99 

1.24 X 10® 

6 

3300880 

161580 

blue 

Dora 

500 

99 

9.92 X 10® 

12 

3997898 

191793 

magenta 

Emma 

1000 

99 

7.94 X lO^o 

24 

4739379 

176633 

red 

Flora 

2000 

99 

6.35 X 10“ 

48 

5046663 

75513 

orange 


Table 1. Features of Le SBARBINE simulations run with PlancklS parameters Qrn = 0.307, r2A = 0.693, ag = 0.829 and h = 0.677 and 
containing 1024^ dark matter particles. The last two columns report the total number of haloes identified with the Spherical Overdensity 
at redshift z = 0 that are resolved with more than 10 and 300 particles, respectively. 


Secondary set of simulations 


name 



crs 

box [h ^Mpc] 

mp[MQh 

colour 

Tea 

0.2 

0.8 

0.7 

150 

1.396 X 10® 

gray-square 

Tea-big 

0.2 

0.8 

0.7 

1000 

4.135 X 10“ 

gray-square 

Tina 

0.2 

0.8 

0.9 

150 

1.396 X 10® 

gray-triangle 

Tina-big 

0.2 

0.8 

0.9 

1000 

4.135 X 10“ 

gray-triangle 

Vera 

0.4 

0.6 

0.7 

150 

2.791 X 10® 

brown-square 

Vera-big 

0.4 

0.6 

0.7 

1000 

8.271 X 10“ 

brown-square 

Viola 

0.4 

0.6 

0.9 

150 

2.791 X 10® 

brown-triangle 

Viola-big 

0.4 

0.6 

0.9 

1000 

8.271 X 10“ 

brown-triangle 

Wanda (wmapT) 

0.272 

0.728 

0.81 

150 

1.898 X 10® 

blue-circle 

Wanda-big (wmap7) 

0.272 

0.728 

0.81 

1000 

5.624 X 10“ 

blue-circle 


Table 2. Details of the small set of 10 simulations with different cosmological parameters. Each contains 512^ dark matter particles 
with initial conditions generated at redshift z = 99. For all the models the Hubble parameter is h = 0.6777, apart from the WMAP7 
cosmology for which h = 0.704. 


But these differences appear only for high u and high red- 
shifts, which play little role in our calibration of the mass 
function parameters. 

All runs were performed in Padova on “Nemo”: a Su- 
perServer Twin 2U Dual Xeon Sandy Bridge composed by 
four independend node servers each equipped with two Xeon 
Sandy Bridge 8 Core E5-2670 and 128 GB of RAM, for a 
total of 64 cores or 128 CPU-threads and 512 GB of RAM. 

2.2 Halo catalogues 

For each stored particle snapshot we identify haloes using 
a Spherical Overdensity (SO) algorithm (e.g. Tormen 1998; 
Tormen et al. 2004; Giocoli et al. 2008). We chose this rather 
than the Friends-of-Friends (FoF) method of Davis et al. 
(1985), because we believe it to be slightly closer to physical 
models of halo formation, and because it is quite similar to 
how mass is defined in observational data. 

For each particle distribution, we estimate the local 
dark matter density at the position of each particle by cal¬ 
culating the distance to the tenth nearest neighbour. In 
this way, we assign to each particle a local density pi oc d~^Q. 
We then sort particles in density and choose as centre of the 
first halo the position of the densest particle. We grow a 
sphere around this centre, and stop when the mean den¬ 
sity within the sphere falls below a desired critical value. At 
this point we assign all particles within the sphere to the 
newly identified halo, and remove them from the global list 
of particles. We then choose the densest particle of the ones 


remaining, and repeat (i.e., grow a sphere around it until 
the mean enclosed density falls below threshold, etc.). We 
continue in this manner until none of the remaining parti¬ 
cles has a local density large enough to be the center of a 10 
particle halo (as we discuss shortly, we apply a more strin¬ 
gent cut when we fit for the mass function); particles not 
assigned to any halo are called ‘field’ or ‘dust’ particles. 

For the critical overdensity we adopt six different defi¬ 
nitions: 2000, 1000, 500 and 200pc(^), 200p6 and the virial 
value. We chose these values of overdensity since they are 
(or they are very close to) the commonly used ones: 200p6 
is motivated by the spherical collapse model in an Einstein- 
de-Sitter universe and - together with 200pc - is a popular 
choice (Tinker et al. 2008); moreover, 200y9c is often used to 
define galaxy cluster masses; 500pc (and the higher overden¬ 
sities) are used in X-ray analyses, and in general in observa¬ 
tions that are able to resolve only the inner parts of haloes. 
The virial overdensity depends on redshift and cosmological 
model (e.g. Bryan & Norman 1998); we use the numerical 
solutions of Eke et al. (1996). 

The comoving density of the background is 

Pcom = Pb — Pc(0)^m(0) — •> (1) 

where Pc(0) = SHq/SttG ~ 2.775 x 10^^h~^MQh^Mpc~^ is 
the critical density. Eor the Planck cosmology adopted in 
the main set of simulations, Avir(z = 0) ~ 319pb is greater 
than 200pb (corresponding to ~ 98pc) and lower than all the 
other thresholds we consider. At high redshifts both 200pb 
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Figure 1. Schematic representation of the halo identifications 
in the particle density distribution at 2 : = 0. Different colours 
represent the various overdensities, from 2000pc down to 200^^. 
Since our SO halo finder starts to grow the sphere starting from 
the densest particle, the halo centered at 1 is common to all the 
catalogues, whereas the ones at 2 and 3 belong only to lOOOpc 
and to the virial catalogue, respectively. This is clearer in the 
right hand panel of the figure which shows the density profiles of 
the three systems. 


and 200y9c converge to the virial definition, the first from 
below and the second from above. 

Fig. 1 shows a schematic representation of the haloes 
identified with our spherical overdensity finder. In contrast 
to previous work, in which an FoF catalogue is used as 
the basis for subsequent SO identifications, we run our halo 
finder code from scratch for each threshold. Hence, although 
a halo in one catalogue may be present in another, this is 
not necessarily true. In particular, while the halo centered 
at 1 is common to all catalogues, the ones at 2 and 3 be¬ 
long only to lOOOpc and to the virial catalogue, respectively. 
This happens, as can be noticed from their density profiles 
presented in the right part of the figure, because while halo 
1 is dense enough to go from overdensity 2000y9c to 200p6 
- so is in common to all catalogues - haloes 2 and 3 reach 
only lOOOpc and the virial overdensity, respectively, and are 
present only in those catalogues. 

Halo identification in each simulation snapshot was 
done with the aim of studying the evolution the halo mass 
function. For each snapshot-density threshold combination, 
we saved a catalogue containing all the information about 
the identified haloes. So as not to be biased by the mass and 
force resolution, we only consider systems resolved with at 
least 300 particles (Maccio et al. 2007, 2008; Velliscig et al. 
2015). Therefore, while the z — ^ catalogues of each sim¬ 
ulation contain many haloes over a wide mass range, only 
the high mass end of the higher redshift outputs is reliably 
measured. 

We also ran an Ellipsoidal Overdensity (EO) finder 
(Despali et al. 2013, 2014; Bonamigo et al. 2015), which we 
used for estimating the triaxial properties of the collapsed 
systems. In what follows we concentrate on the results for 
SO haloes; a brief summary of the corresponding results 
for ellipsoidal haloes can be found in Appendix A. These 
are broadly similar to spherical haloes, although the best fit 
mass function parameters differ slightly from those for the 
SO case. 


3 MODEL FOR THE HALO MASS FUNCTION 


Let dn/dlnM denote the comoving number density of 
haloes in a logarithmic bin dlnM around mass M. Then, 
the mass fraction in such haloes is 

/(M)dM = —-^ArdlnM. (2) 

^ Pf, dlnM ^ ^ 

It is usual to define 

(7^ (M, z)=J^ k^Pnnik,z) ^2 , (3) 

where Piin{k,z) is the initial power spectrum extrapolated 
to redshift z using linear theory, W{x) = 3ji{x)/x, and 
R{M) is given by requiring M/(47rR^/3) = pb (recall that 
Pb is comoving, so it is independent of redshift). In principle, 
a{M,z) depends only on P\in{k,z) and the smoothing filter 
W. In practice, the number of Fourier modes that can be 
effectively sampled in the initial conditions depends on some 
computational limits (e.g. box size, number of particles, etc). 
As a result, there are differences between the actual power 
spectra in a box and the theoretical mean value. Appendix 
B illustrates the impact on 

We compute for each box using the actual re¬ 

alisation of P(k) in it. This reduces bias and scatter in the 
mass function, especially in the high mass tail of each simu¬ 
lation, allowing us to reach higher precision (down from ten 
to a few percent) particularly in the smaller boxes within 
which cosmic variance would otherwise contribute substan¬ 
tially. In particular, the difference between the theoretical 
input power spectrum and the one measured from each sim¬ 
ulation can be seen in Figure Bl, while the associated dif¬ 
ference in the mass function is shown in the bottom panel 
of Figure B2. Appendix B also discusses the impact of box 
size on the mass functions, and compares a number of differ¬ 
ent prescriptions that account for finite box-size effects on 

Although a number of workers so far have parametrised 
the mass function in terms of a alone (i.e., they write 
/(M) dM = /(cr)dcr and work with /(cr)), Sheth & Tormen 
(1999) were careful to parametrize in terms of 

V = 5l{z)la^{M), (4) 

where Sc{z) is the critical linear theory overdensity Sun (z) 
required for spherical collapse divided by the growth factor 
(Carroll et al. 1992). This quantity, which depends weakly 
on Q for the ACDM family of models, is well-approximated 
by 


Suniz) « T (i2n)P^[l + 0.0123 Ign(z)] (5) 

(Kitayama & Suto 1996). The rationale for including it came 
from the fact that they used the virial density to find haloes, 
and the same model which predicts this virial density also 
predicts Sc- Court in et al. (2011) highlight the fact that this 
factor becomes increasingly important at high masses; fail¬ 
ure to include it can masquerade as non-universality. 

For this reason, our reference parametrisation for the 
halo mass function is that of Sheth & Tormen (1999): 

f{M)dM = f{iy)diy (6) 


where 



1 

ly'p 



e 


—1/ i‘i 




( 7 ) 
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Figure 2. Halo mass function at four redshifts, for all SO haloes with more than 300 dark matter particles in the SO catalogues, 
identified using the redshift dependent virial density. Green, cyan, blue, magenta red and orange symbols show results for Ada, Bice, 
Cloe, Dora, Emma, and Flora. Black dashed line, same in all panels, shows the result of fitting the z = 0 points to Eq. (7). The middle 
panels show the residuals (in log space) from this best fit. The dashed vertical lines show the maximum i/ for which the bins contain at 
least 100 objects. The two solid vertical lines show the minimum u at which we may expect differences of at least 5% in the halo mass 
function accounting for 2LPT initial condition, adopting the rescaling of equation (11) by Reed et al. (2013), for the simulations starting 
at 2 ; = 99 and 2 ; = 124 (right and left lines, respectively). Red solid line shows the best fit mass function obtained using the data from 
all the snapshots with 2 : ^ 1.25. The bottom panels show residuals from the fit to the scaled counts from all snapshots up to 2 : = 1.25. 
The departures at high redshift are reduced, so universality is even more pronounced. 


with u = aiy. The parameters (a,p, Aq) define the high-mass 
cutoff, the shape at lower masses, and the normalisation of 
the curve, respectively. In addition, because P\in{k) deter¬ 
mines the value of the mass variance S{M) = which 

enters in the definition of ly, the initial power spectrum plays 
an important role. Since it enters in the denominator of ly, 
one might say that the mass function is ‘non-perturbative’ 
in Piin{k). 


4 THE UNIQUENESS OF THE VIRIAL 

OVERDENSITY 

This section studies the universality of the halo mass func¬ 
tion for haloes identified using the virial overdensity. We do 
so by finding the set of parameters (a,p,Ao) that best fit 
the z = 0 data from the PlancklS simulations. We then use 
measurements at other redshifts and cosmologies to test for 
universality. We start with the virial overdensity because 
we believe it to be the most physically motivated choice for 
identifying haloes. Later in the paper we study haloes iden¬ 
tified using the other SO catalogues. 

Our model for the mass function, equation (7), has 
three free parameters (a,p, Ao), whose values we adjust so 
as to minimise the Chi square with respect to the measured 


binned mass function: ^ 

X^{a,p,Ao) = ^ - l°g(^/(^))/»«)' (8) 

where the sum is over binned counts. The bins are equally 
spaced in log 3 ^Q(z/), with Alog 3 ^Q(z/) = 0.05 and we neglect 
covariances between the binned counts when fitting. ^ We 
discard bins with fewer than 30 objects (typically the high- 
ly bins), to set a limit on the Poissonian error eiog(iy/(iy)) in 
a bin. There were no significant differences in the best-fit 
parameters when we repeated the analysis using only bins 
with at least 100 objects. 


^ Note that Sheth &: Tormen (1999) only varied a and p: they 
fixed the value of A by requiring that the integral over all masses 
give pfo. 

^ We do not include a bin that counts the mass fraction which is 
not assigned to haloes. See Manera et al. (2010) for an algorithm 
which does not require binned counts, and does account for the 
unassigned mass fraction. 
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4.1 The halo mass function at different redshifts 
and overdensities 

4.1.1 Virial haloes 

The best-fit values for the z = 0 virial overdensity halo 
counts are 

{a,p,Ao) = (0.794±0.005,0.247±0.009,0.333±0.001). (9) 

The dashed black curves, same in all panels of Figure 2, show 
equation (7) with these parameters. The symbols in the dif¬ 
ferent panels show the measured virial halo mass function 
at four different redshifts as labelled. Green, cyan, blue, ma¬ 
genta, red and orange points show results from the six sim¬ 
ulations of the main set: Ada, Bice, Cloe, Dora, Emma and 
Flora. Note the excellent agreement among the simulations 
despite the very different volume that each samples. This is 
because the way we calculate iy{Sc, a) largely eliminates the 
impact of cosmic variance. 

The middle panels of Fig. 2 show logarithmic residuals 
from the best-fit to the z = 0 counts. These indicate the 
goodness of our fit, which was calibrated in log space. The 
residuals are close to zero at small and intermediate values 
of 12 . They are larger at high z/, where the best-fit predicts 
more haloes than we find in the simulations at z = 0. Some of 
this is due to the Poisson uncertainty in the small number 
counts; at higher redshifts the effect is reduced since the 
high 12 tail becomes populated by lower mass haloes. The 
smallness of the residuals in the other panels indicates that 
the z = 0 model is a good description of the virial halo mass 
function at higher redshifts as well. 

The vertical lines in the bottom panels show the regime 
of influence of two numerical effects: (i) First, the dashed 
gray vertical line shows the maximum 12 for which there are 
at least 100 (rather than 30) haloes in the bin. This mat¬ 
ters only at high redshifts and in the high mass regime and, 
as mentioned before, the best fit parameters are not sig¬ 
nificantly different; (ii) Second, we addressed the effect of 
not using 2LPT initial conditions for our simulations: we 
rescaled our points using equation (11) of Reed et al. (2013) 
(calibrated only for 1 ^ ^ 5), which gives an estimate of 

the bias in the halo mass function between ZA and 2LPT 
methods; the two solid vertical lines show the minimum i/ at 
which we may expect differences of at least 5%, for the sim¬ 
ulations starting at z = 99 (right) and z = 124 (left). Such 
discrepancies would mainly affect the high-; 2 ; and high-i/ data 
points; they do not affect the counts in the z ^ 1.25 range we 
will adopt to calibrate the parameters of equation (7). The 
rescaling equation by Reed et al. (2013) and the ZA/2LPT 
difference are discussed in detail in Appendix B4. 

From the figure and the discussion above we conclude 
that there is no significant systematic deviation from uni¬ 
versality within 8% for 12 < 10 at redshifts z ^ 5. Therefore, 
we can combine measurements at many different redshifts 
to increase the precision in the estimate of (a,p, Ao). In ad¬ 
dition to increasing the statistics, adding the high-redshift 
data fills-in the h.igh -12 tail, allowing it to play a greater 
role in determining the best-fit parameters. We determined 
best-fit parameters using the virial halo data from all the 
snapshots between z = 0 and 1.25 (a total of 15 snapshots). 
At ^ 1.25, the data no longer samples the whole mass 
function, preventing a reliable fit. The resulting best fit pa- 
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Figure 3. Change in the best-fit mass function when one param¬ 
eter at a time is changed by up to ±30% (in increments of 10%). 
Blue and red curves show cases where the parameter values are 
larger or smaller than the reference one. 

rameters are: 
a = 0.7663 ± 0.0013 

p = 0.2579 ± 0.0026 (10) 

Ao = 0.3298 ± 0.0003. 

The differences between this all-z fit and the previous z = 0 
fit are of order Aa ~ 3%, Ap ^ 4% and AAo ^1%. This 
leads to percent level differences in the mass function which 
reduce the residuals at high redshift. The red solid line in 
each panel of Figure 2 shows this all-z best fit: it can hardly 
be distinguished from that for z = 0 only (black dashed 
curve) but, as shown by the residuals in the lower panels, it 
traces the high -12 part of the mass function more closely. 

Figure 3 shows how the virial mass function changes as 
each of the parameters (a,p, Aq) is varied by ±30% while the 
other two are kept fixed. While the mass function is sensitive 
to both Ao and a - which modify the normalisation and the 
high mass cut-off, respectively - this is less true for p. Even 
if p changes by ±30%, the mass function changes by 5% at 
most, for 12 < 0.1 and 12 > 10. In what follows, we use the 
values in equation (10) as our reference model. 

4.1.2 Other overdensity thresholds 

We now study the mass functions associated with the other 
five density thresholds. It is tempting to identify these 
thresholds with an ‘effective formation redshift’, and hence 
with an effective value of 6c(z) when defining 12 . Since the 
role of a in equation (7) is simply to rescale u, one might 
wonder if these other halo definitions lead to universal mass 
functions which differ from those for the virial overdensity 
only in the value of a. 

Table 3 shows the values of the best-fitting parameters 
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at z = 0 for the different threshold densities. Clearly the 
other parameters, p and Aq, also depend strongly on how 
the haloes were identified. Moreover, as we will show below, 
universality clearly does not hold for any of these other def¬ 
initions: in all cases, the best fit parameters for the higher 
redshift counts depart significantly from those at z = 0. 

To illustrate this. Figure 4 shows the halo mass func¬ 
tions at z = 0, 1, 2, and 5 for haloes identified using 200pb 
as density threshold. The dashed curve, which is the same 
in all four top panels, shows the best-fit to the z = 0 data 
points. The middle panels, which show the residuals with re¬ 
spect to this fit, show that it overestimates the counts at all 
higher redshifts. These departures from universality are in 
agreement with previous work on 200pb haloes (e.g. Tinker 
et al. 2008). The solid curves in the upper panels show the 
result of rescaling our universal virial counts as described in 
the next section; the bottom panels show the residuals with 
respect to it, which provides a much better fit. 

Figure 5 shows a similar analysis of haloes identified 
using 200pc(^)- In this case, the best-fit relation at z = 0 
underestimates the counts at higher z. Finally, Fig. 6 shows 
haloes identified using threshold values of 500, 1000, and 
2000pc(^), at redshifts z = 0 and z = 1. The trends are 
qualitatively similar to those for 200pc(^), with the z = 0 fit 
underestimating the counts at higher z. 


4.2 The halo mass function for different 
cosmologies 

Having established that the virial mass function is universal 
with respect to redshift, we now test if it is universal across 
other background cosmological models. To do so, we use the 
halo catalogs in our secondary set of 512^ simulations, whose 
properties are summarised in Table 2. 

First, for virial haloes, we used the best fit parameters 
of equation (10) which were calibrated using the 1024^ par¬ 
ticle simulations of a Planck 13 cosmology. The four panels in 
Figure 7 show the logarithmic difference from this best fit at 
redshifts z = 0,1, 2, 5. Following the colour code of Table 2, 
gray points are for Qrn = 0.2, brown for Qrn = 0.4 and light 
blue for Qrn = 0.272 (the WMAP7 cosmology); squares and 
triangles represent ag = 0.7 and ag = 0.9, respectively. For 
halo mass definitions that are close to the virial value the 
residuals are relatively small, indicating that our relation - 
calculated exclusively from the Planck 13 cosmology data - 
remains valid for the other cosmological models as well. In 
addition, although we do not show this explicitly, for these 
cosmological models also, universality is broken when con¬ 
sidering other overdensities. This demonstrates universality 
of the virial relation with respect to other cosmological mod¬ 
els. However, departures from this universality may arise if 
one considers more extreme departures from Planckl3. E.g., 
(Courtin et al. 2011) explore models in which the growth fac¬ 
tor, growth history and virial density differ more radically 
from that of the Planck cosmology, and find correspondingly 
larger departures from universality. 

Figure 8 shows the results from all our simulations and 
redshifts together. Different colours represent the different 
simulations (with the same colour code as in previous fig¬ 
ures), and different symbols distinguish the four redshifts 
(and not the different simulations). The best fit to all the 





Figure 6. Same as Fig. 2, but for haloes identified using 500, 
1000 and 2000 pc instead of at 2 ; = 0 and 1. 
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p(SO) 

a 

P 

A 

z = 0 

200pb 

0.739 ±0.005 

0.206 ±0.008 

0.360 ±0.001 

^vir 

0.794 ±0.005 

0.247 ±0.009 

0.333 ±0.001 

200pc 

0.903 ±0.006 

0.322 ±0.009 

0.287 ±0.001 

500pc 

1.166 ±0.009 

0.344 ±0.012 

0.236 ±0.001 

lOOOpc 

1.462 ±0.012 

0.349 ±0.015 

0.197 ±0.001 

2000pc 

1.821 ±0.017 

0.413 ±0.017 

0.158 ±0.001 

All z - Planck cosmology 

^vir 

0.7663 ±0.0013 

0.2579 ±0.0026 

0.3298 ±0.0003 

All z & cosmologies 

^vir 

0.7689 ±0.0011 

0.2536 ±0.0026 

0.3295 ±0.0003 


All z &: cosmologies - Cluster Counts: M. 

yir > 3 X IO^^Mq/Zi 

^vir 

0.8199 ±0.0010 

0 

0.3141 ±0.0006 


Table 3. Dependence of best-fit parameters on the overdensity used to identify SO haloes. In the top part we show the three parameters 
calculated at 2 ; = 0; then we report those obtained by fitting the virial halo counts from all snapshots up to 2 ; = 1.25 {i) of the Planck 
simulations and {ii) of all the simulations with different cosmologies together. The bottom row shows the parameters which best-fit the 
^vir ^ 3 X 10^^ counts from all cosmologies and all redshifts. 



1/ u 1/ 1/ 


Figure 4. Same as Fig. 2, but for haloes identified using 200p5 instead of Ayir. Middle panels show residuals with respect to the best-fit 
at 2 : = 0 (dashed curve, same in all top panels); these show that the mass function is not universal across all redshifts. Lower panels show 
residuals with respect to our rescaled model (equation (12); solid curves, different in each top panel). 


virial halo catalogues out to z = 1.25 has 
a = 0.7689 zb 0.0011 

p = 0.2536 zb 0.0026 (11) 

^0 = 0.3295 zb 0.0003, 

and is shown by the solid black line. Lower panel shows 
the logarithmic residuals. These best-fit parameter values 
equal - at the per mil level - the values obtained in the last 


Section from our first set of simulations, which assumed a 
PlancklS cosmology (equation 10), thus further confirming 
the universality as a function of the background cosmological 
model. 

Table 3 compares the three sets of best-fitting values 
for the virial haloes. The bottom row shows the parameters 
which best fit the counts of cluster-mass haloes: Myir ^ 
3 X 10^^h~^ Mq. Notice that at these high masses p = 0, 
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Figure 5. Same as Fig. 2, but for haloes identified using 200pc; middle panels show that the universality is broken in the opposite sense 
to when the threshold was 200^;,, and bottom panels show that our rescaled model accounts for this quite well. 


z = 0 


z = 1 


z = 2 


z = 5 



Figure 7. Cosmology independence of the virial halo mass function. We show the logarithmic difference from the best fit calculated 
using the Planck cosmology for redshifts 2 ; = 0,1, 2, and 5. Following the colour code of Table 2, gray points indicate the results of the 
simulations with Qrn = 0.2, brown for Qrn = 0.4 and light blue for Qrn = 0.272; squares and triangles represent ag = 0.7 and ag = 0.9, 
respectively. 


SO that the mass function is like that of (Press & Schechter 
1974). We discuss this further in Appendix C in the context 
of cosmological constraints in the flm-crs plane from cluster 
counts. 

Covariances between the best-fit parameters are shown 
in Fig. 9. Contours show 1-, 2- and 3-a levels for each pair 
of parameters. While Aq and a are not strongly degener¬ 
ate, p correlates with both a and Aq. Manera et al. (2010) 
show that such correlations can be understood as resulting 
from requiring the model to reproduce the total measured 
mass fraction in haloes (the mass fraction is measured with 
much greater precision than is the detailed shape of the mass 
function; also see Figure 3). 


5 RESCALING THE MASS FUNCTION: A 

UNIVERSAL PARAMETRISATION 

We now describe a simple method that allows one to derive 
the halo abundances associated with other overdensities at 
any redshift by a straightforward rescaling of the virial halo 
mass function. 

We calculated the best-fitting parameters for all the six 
halo catalogues at all outputs between z = 0 and 1.25. For 
the overdensities ^ 500pc(^) we only fitted the data points 
out to z ~ 0.4 so as to be sure of having good statistics 
for both the shape and the cut-off of the mass function, i.e. 
p and a, respectively. At higher redshifts the small number 
of haloes which probe a smaller range of u cannot break 
degeneracies in the determination of the three parameters. 
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Figure 8. Universality of the virial halo mass function. Symbols show measurements from all our simulations - both from the main and 
secondary set - at four redshifts: z = 0,1, 2, 5. Although we show measurements at these four redshifts, the fit, which was calibrated on 
all simulations, used only the 2 : ^ 1.25 snapshots. 





Figure 9. Covariance between fitted parameters of the virial halo mass function (the 1024^ and 512^ runs, and all redshifts to 2 : ^ 1.25): 
Contours show 1-, 2- and S-cr reference levels for each pair of parameters. 


Nevertheless, even after restricting to z ~ 0.4, we obtain 
very good results. As the reference virial halo mass function 
we used the fit from eq. (11). Fig. 10 shows how the best fit 
parameters vary as a function oi x — log(A( 2 :)/A-f;ir (4). 

As noted in Tinker et al. (2008) (albeit for a differ¬ 
ent functional form), the best fit parameters (a,p, Aq) are 


smooth functions of the critical density threshold. As can 
be seen in the top panel of Fig. 10, the normalisation Aq 
decreases linearly with the overdensity - a natural conse¬ 
quence of the decrease in halo mass - while a and p both 
increase with threshold overdensity. The trends can be de- 
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scribed using linear or quadratic polynomials: 

a = 0.4332 + 0.2263 a: + 0.7665, 

p = -0.1151 a:^ + 0.2554a:+ 0.2488, (12) 

A = -0.1362 a:+ 0.3292, 

shown by the blue solid lines in Figure 10. While a and A 
behave very regularly, the determination of p is less certain. 
This does not have a big influence on the rescaling method, 
since the mass function is less sensitive to p than to the 
other parameters (see Figure 3). The red dashed curves in 
Figure 10 show the trends obtained by fitting to the z = 0 
counts only: 

a = 0.3881 xl + 0.2776 xo + 0.7837, 

p = -0.07459a:o + 0.2016a:o +0.2518, (13) 

A = -0.1337 a:o+ 0.3315, 

where xo = log(A(zo)/Ayir(zo = 0)). These trends are very 
similar to the previous ones. 

The above relations show that both a and p in¬ 
crease with increasing threshold, qualitatively consistent 
with Manera et al. (2010) who studied FOF rather than 
SO haloes. (They found that a and p increase as the FOF 
linking length is decreased, and it is well known that shorter 
linking lengths return denser haloes.) In addition, upon not¬ 
ing that a multiplies Sc{z) in Eq. (7), the increase of a with 
density threshold is qualitatively consistent with the notion 
that the denser inner parts of a halo virialized at higher 
redshift, so mass functions for higher density thresholds are 
qualitatively like those for higher redshifts. Unfortunately, 
the agreement is not quantitative. To see this, consider the 
thresholds 2000pc and 200pc which differ by a factor of 10. 
If the mass associated with the denser threshold virialized 
at (1 + Z 2 ooo)^ ~ 10, we would expect the associated Sc to 
be larger by a factor of (1 + ^ 2000 ) ~ 10^^^. If this increase 
is to be provided by increasing a, then a must be larger by 
102/3 jg considerably larger than the ratio of a values 

for 2000pc and 200pc in Table 3. 

Even though we do not have a quantitative physical 
understanding of the scaling in equations (13) and (12) - 
a question we return to in Section 6 - they can still be 
used to predict the mass functions at overdensities of current 
interest that are not considered in our work. These precise 
fitting formulae are able to follow the evolution of the non- 
virial mass functions at all redshifts, at the level of a few 
percent - which is comparable to the intrinsic uncertainty in 
the mass function. 

To show this, we return to the solid lines and the bot¬ 
tom panels in Eigure 4-6. The black solid curves in Eig. 4-6 
show equation (12) (recall that the dashed curves show the 
best z = 0 fits) and the bottom panels show the logarithmic 
difference from our model. The agreement is indeed very 
good: the residuals are generally smaller than a few per¬ 
cent; this is quite acceptable given the intrinsic scatter of 
the mass function is of this order. The middle and bottom 
panels show that our model curves are able to account both 
for the change in the normalisation and the tilt of the mass 
function as z changes. 

To conclude: our results provide an efficient way for pre¬ 
dicting the halo mass function for any redshift and overden¬ 
sity in a Planck ACDM-like model. Once log(A(z)/A-i;iT. (z)) 





Figure 10. Dependance of best-fit parameters a^p^Ao on x = 
\og(A(z)/Avir{z)) for 2 ; = 0 to 1.25. Different coloured sym¬ 
bols represent different overdensities; smooth curve shows equa¬ 
tion (12). 


has been calculated, it is straightforward to obtain the other 
parameters. 

Returning to the other cosmological models, Eigure 11 
shows the logarithmic differences from our reference rescal¬ 
ing models for all the mass functions of the secondary set, 
at redshifts z = 0,1. At both redshifts, the residuals are 
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Figure 11. Mass functions for other cosmological models: Loga¬ 
rithmic difference from our best fit models in the secondary set of 
512^ simulations, at two different redshifts {z = 0,1) and for 
all the overdensities. Gray, brown, light blue show results for 
= 0.2, 0.4 and 0.272. Squares and triangles show results for 
(78 = 0.7 and 0.9. For the virial haloes we use the fit obtained 
from stacking all 2 ; outputs of the Planck cosmology. For the other 
cases, the residuals are with respect to equation (12). 


very small showing that our fitting relations also work for 
any cosmological model once mass and redshift are appro¬ 
priately rescaled to v. 


6 MATCHED HALOES 

So far, we used the halo catalogues derived independently 
at each overdensity threshold from the particle density dis¬ 
tribution saved in the snapshot files. As already noted, the 
number of haloes identified using different overdensities at a 
certain redshift is not the same. This means that we are not 
looking at a “rescaled version” of the virial population, but 
at different halo samples (see Figure 1). For example, virial 
haloes may contain many smaller higher density peaks, each 
of which could correspond to a massive substructure within 
the virial radius. On the other hand, the 200p6 threshold is 
lower than the virial one and so - at low redshifts - the iden¬ 
tified haloes will be larger and more massive then the virial 
ones. In general in this latter case we may have fewer haloes 
than in the virial catalogues, since the particles correspond¬ 
ing to more than one smaller virial halo may be included in 
a single 200p6 halo (c.f. the virial halo 3 in Figure 1). 

In this section we analyse the halo mass function of 
“matched haloes” which we create as follows. For each virial 
halo, we select from the different halo overdensity catalogues 
the object whose center of mass is closest to that of the virial 
halo. For the halos identified with a threshold that is denser 
than virial, we only keep the object which corresponds most 
closely to the virial one. To ensure good resolution (and 
hence a match) in the inner regions, we only considered virial 
haloes with more than 10^ particles. In this way, we define 
a series of approximately concentric spheres of ever higher 
density within the virial radius: a density profile. 

If we average together all the profiles for a narrow bin 
in virial mass, then we can use this mean profile, and the 
scatter around it, to model the ‘matched halo’ mass function 
associated with different overdensity thresholds A ^ 
because 

If p{M/:^\Mvir) is sharply peaked around 
Q/\ (A/7)7.-r) say, then 

d?7//\ ^Tlyi'p dTVL-yjT^ 

dMA dMvir (IMa 

so that 

as in Bocquet et al. (e.g. 2015). But if the scatter around 
the mean profile is not negligible, then the full convolution 
of equation (14) must be performed. 

In principle, p{M/:^\Mvir) could be determined directly 
by fitting an NFW profile Navarro et al. (1996) to each 
virial halo in the stack. This functional form has just one 
free parameter, the concentration c, so p{M/:^\Myir) is sim¬ 
ply related to the mean c at fixed Myir^ and the scatter 
around this mean, which is known to be lognormal with 
width (7inc|m = 0.25. 

In practice, we have not measured the profile shapes. 
Rather, we use the fact that the concentration of a halo is 
related to its mass accretion history (Zhao et al. 2009), which 


(14) 
Ma = 

(15) 

(16) 
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we also do not measure. Rather, we estimate it, and hence 
c, using the model of Giocoli et ah (2012c). We then use the 
mean c along with the assumption that the scatter is negli¬ 
gible, to compute the value of qa to insert in equation (16) 
for the desired A. 

The symbols in the left panels of Fig. 12 show the mass 
function of the haloes in the different catalogues matching 
the virial systems, at redshift z — 0 and z = 1 for the cat¬ 
alogues constructed considering 200, 500, 1000 and 2000pc 
as thresholds. The solid curves, following the corresponding 
data points, show the virial mass function rescaled to the 
different overdensities. 

The panels on the right of Figure 12 show the logarith¬ 
mic residuals of the data points with respect to the corre¬ 
sponding rescaled models. They all remain well below 10% 
- apart from the high-i/ tail - indicating that the preci¬ 
sion of the rescaled mass functions is of the same order as 
of the “true” virial mass function. Some of the discrepancy 
may be caused by the fact that not all haloes follow NFW 
profiles all the way down to the very center (Einasto 1965; 
Retana-Montenegro et al. 2012; Ludlow et al. 2013; Dutton 
k Maccio 2014). 

The case of 200pb must be treated separately, since the 
more massive haloes identified using 200pb may include more 
than one smaller virial overdensity haloes in their outskirts 
which we exclude (e.g. object 3 in Figure 1). Nevertheless, in 
this case also, the rescaled relation captures the behaviour 
of the mass function measured in the simulations with a 
precision comparable to the other overdensities. Recall that 
the rescaled mass function at 200pb at z = 1 is almost the 
same as the virial mass function, since Avir at high redshift 
is nearly equal to 200p6. 

Figure 13 shows how different these rescaled mass func¬ 
tions are from the original fits of Table 3. The upper panels 
show the two mass functions - the rescaled one in dashed 
lines and the original ones in solid lines - for all the den¬ 
sity thresholds; the lower panel shows the differences in halo 
counts between the two cases. Notice that the matched mass 
functions ’lose’ haloes: because some of the smaller denser 
systems are halo substructures. This is true for all the over¬ 
densities except for 200pb: as we said, in this case we are 
outside the virial haloes and so we may in fact have more 
haloes than those found at the virial overdensity. 


7 COMPARISON WITH PREVIOUS WORK 

Figure 14 shows the residuals between our virial mass func¬ 
tion lyfiiy) - the one calibrated from all z ^ 1.25 outputs of 
the 1024^ and the 512^ simulations - and the mass functions 
derived in previous work (as indicated). The shaded region 
in each panel shows the range in z/, at z = 0, over which the 
various authors have calibrated their mass functions. The 
two vertical dot-dashed lines show the range in ly over which 
we have performed our calibration. 

To facilitate comparison, the main features of these 
works are listed in Table 4. Most of these authors used 
a FOF algorithm to identify the haloes in their simula¬ 
tions: while b=0.2 is the most common choice for the link¬ 
ing length, Manera et al. (2010) use three different values 
(b=0.15, 0.168, 0.2); moreover, Warren et al. (2006) discuss 
how to correct the FOF masses, claiming that the identifi- 
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Figure 13. Difference between the mass functions of Table 3 
(solid lines) and the rescaled ones from the matched haloes as 
in Figure 12 (dashed lines); the lower panel shows the residuals 
between the two relations on the halo counts 


cation is influenced by how well an halo is resolved and thus 
the same linking length may identify structures with differ¬ 
ent enclosed overdensities. The same argument is discussed 
also in Courtin et al. (2011) and More et al. (2011) and they 
predict a linking length equivalent to our virial overdensity 
(at z=0) of 0.193 and 0.206, respectively. Sheth k Tormen 
(1999) and Tinker et al. (2008) use SO algorithms, while 
Watson et al. (2013) compares the results of the two meth¬ 
ods (in Figure 14 we show only their universal FOF fit). At 
intermediate masses all the analytical mass functions agree 
quantitatively. Larger differences arise for more massive sys¬ 
tems: in this range, the precision of the fit is strongly affected 
by the resolution of the simulation and consequently by the 
number of high mass haloes. 

The main purpose of our work is to analyse the impact 
of the density threshold chosen to identify the haloes on the 
universality of the mass function: the majority of the iden¬ 
tification algorithms use a threshold A of 178 or 200 for all 
redshifts (e.g. Tinker et al. (2008) and Watson et al. (2013)). 
Since we have argued that ignoring the redshift evolution of 
A causes most of the observed non-universality one must 
carefully match density thresholds when comparing differ¬ 
ent works. This is shown more clearly in the lower panels of 
Figure 14, where we compare our results with those of Tin¬ 
ker et al. (2008): we tried to match their SO thresholds as 
closely as possible to ours, to test their compatibility. The 
Figure shows - from top to bottom - the residuals between: 
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Figure 12. Top — left: z = halo mass functions of the matched haloes. Different colours show different overdensities, while the black 
solid curves show the corresponding rescaled mass functions, obtained from Eq. 16. Top — right: logarithmic residuals from the rescaled 
mass functions of the upper panel. The scatter with respect to the rescaled mass function is of the same order of the intrinsic scatter of 
the mass function (shown by the virial case). Bottom: Same as top, but for z = 1. 


(i) their A = 200 fits with our virial fits; (ii) their A = 300 
fits with our virial fits; (Hi) their A = 200 fits with our 
200pb fits. The agreement for cases ii and in is much bet¬ 
ter than for i. Since our virial overdensity at z=0 is closer 
to A = 300 than A = 200, this shows that different anal¬ 
yses converge provided one compares compatible halo iden¬ 
tification schemes (Knebe et al. 2013). Moreover, not only 
are our 200pb fits close to the main result of Tinker et al. 
(2008), but the non-universality they report is compatible 
with the redshift evolution of our 200pb mass function (Fig¬ 
ure 4): the systematic deviation from universality is primar¬ 


ily due to the choice of the overdensity threshold used in the 
halo identification method. The agreement, once overdensi¬ 
ties have been correctly matched, is even more reassuring 
when one notes that Tinker et al. (2008) fit their counts to 
a functional form that differs from our equation (7). Addi¬ 
tional reassurance comes from the fact that, while one may 
expect trends arising from differences in halo identification 
methods, our SO results lie in between the I = 0.168 and 
I = 0.15 FOF results of Manera et al. (2010), just as they 
should. 


As a final remark we note that, even when the compu- 
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tational schemes are theoretically similar, the mass function 
for structures identified with different codes (in the same 
simulation!) can differ by up to 10% (Knebe et ah 2011). 
This might explain some of the discrepancies between our 
best fit model and those of other works. Finally, we underline 
that there could be other effects depending on the assumed 
cosmological parameters as explored by Murray et al. (2013). 


8 DISCUSSION AND CONCLUSIONS 

In this paper we analysed a set of cosmological simulations 
in order to study the halo mass function and its dependence 
on redshift, cosmology and the halo identification method. 
In what follows we summarise our main results. 

(i) We demonstrated and confirmed the universality of 
the virial halo mass function, by comparing the measured 
mass function at many different redshifts and in simula¬ 
tions with different cosmological models. As stated in Sheth 
& Tormen (1999); Courtin et al. (2011), the halo mass 
function at any time, and for any cosmological model, can 
be well described by a single functional form once mass 
and redshift are appropriately parametrised in terms of 

= 5l{z)/c7\M). 

(ii) We showed that only the virial overdensity leads to 
a universal halo mass function: most of the non-universality 
seen in other works arises from not using the virial value 
when identifying halos. Commonly used values of 178 
or 200 X the critical or background density induce non- 
universal trends. 

(iii) We derived three different sets of best fit parameters 
for the virial halo mass function: 

• a — 0.794, p = 0.247 and A = 0.333 - from the z=0 
virial mass function in the Planck cosmology; 

• a — 0.7663, p = 0.2579 and A = 0.3298 - from the 
virial mass function in the Planck cosmology, using all the 
data points up to z = 1.25 (15 snapshots x 6 runs); 

• a — 0.7689, p = 0.2536 and A = 0.3295 - from the 
mass functions in all cosmologies, using all the data points 
up to z = 1.25 (15 snapshots x 16 runs). 

The last two are in remarkable agreement - at the per mil 
level for the best-fit parameters, and sub percent for the 
mass function), illustrating the level of the universality as a 
function of cosmology. In general, the normalisation is the 
most stable parameter, while a and p change between the 
first and the other two cases: when using only the z = 0 
points, the high-z/ tail is less well resolved, and so the fit 
is less precise in the determination of the full shape of the 
mass function. 

(iv) We presented a simple rescaling method which allows 
one to estimate the three parameters of the fitting function 
through first or second order scaling relations. The three pa¬ 
rameters {a,p,Ao) are smooth functions of the overdensity, 
for any redshifts. Equations (12) are able to describe the 
change in slope and normalisation of the mass function (as 
a function of redshift and overdensity) with good precision 
(Figs. 4, 5, 6 and 11). 

(v) Finally, we studied the mass function of “matched 
haloes”: the counterparts of virial haloes at different over¬ 
densities. We provide an efficient rescaling method, equa¬ 
tion (16), which uses knowledge of the mass density profile 




Figure 14. Comparison with some previous work. In all panels, 
the dotted vertical gray lines show our range in mass (for virial 
haloes at 2 ; < 1.25) and the gray region show an estimate of the 
range in mass of the other works; we calculated this last convert¬ 
ing their mass ranges in i/, so it may not be perfectly the same but 
it gives an estimate of the overlapping regions. Top: we show the 
relative residuals of our virial mass function with respect to 

the models of other works. From top to bottom, we compare our 
model with Sheth & Tormen (1999), Jenkins et al. (2001), War¬ 
ren et al. (2006), Watson et al. (2013) and Manera et al. (2010). 
For this last, we show the results obtained from FOF catalogues 
with linking lengths of 0.2, 0.168 and 0.15, which they fit to the 
same functional form we do. Our results lie between their fits to 
the smaller linking catalogs, as they should, over the whole range 
of masses. Bottom: we compare our results with those of Tinker 
et al. (2008) at different overdensities. First we show the relative 
residuals of their A = 200 reference model with our best fit cal¬ 
culated at the virial density; secondly we compare their model 
with our best fit at 200pi). Finally, we test the opposite situation, 
comparing our virial model with their mass function at A = 300, 
an overdensity closer to the value of Ayir adopted in this work. 
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Popart 

box [h ^ Mpc] 

cosmological models 

algorithm 

threshold 

Sheth & Tormen 1999 

(3 sim x) 256^ 

85, 141 

S/O/ACDM 

SO 

y 

< 

II 

<1 

Jenkins et al. 2001 

256^, 512^, 10^ 

84 -3000 

t/ACDM 

FOF 

b=0.2 

Warren et al. 2006 

10243 

96 - 3072 

ACDM 

FOF 

b corrected 

Watson et al. 2013 

30723-60003 

11 - 6000 

ACDM (WMAP5) 

FOF (SO) 

b=0.2 (A=178) 

Manera et al. 2010 

(49 sim x) 640® 

1280 

ACDM 

FOF 

1=0.15, 0.168, 0.2 

Tinker et al. 2008 

5123-10243 

80 - 1280 

ACDM (WMAPl-3) 

SO 

A=n X 


Table 4. Schematic (i.e. incomplete) summary of the main features of the other works to which we compare our results in Figure 14 
(Sheth &: Tormen 1999; Jenkins et al. 2001; Warren et al. 2006; Watson et al. 2013; Manera et al. 2010; Tinker et al. 2008). We list the 
resolution and scale of their simulations and the main algorithm with the associated threshold used to identify the haloes. 


and the concentration-virial mass relation, to estimate their 
mass function with a high precision (Fig. 12). 

These rescaling methods are useful for comparing analyses 
which use different definitions of halo mass. In particular, 
they can be used to translate our universal virial halo mass 
functions to the nonuniversal form associated with halo def¬ 
initions which are closer to those commonly used in observa¬ 
tional studies - ranging from the X-ray and SZ to the visible 
and near infrared. 

We conclude that - over the range of redshifts and cos¬ 
mological models in our simulation set - the virial halo mass 
function is, to within 5 — 8% for a wide range of masses, a 
universal function of redshift and cosmology. Non-universal 
behaviour can be an artifact induced by the halo identifica¬ 
tion method and by the choice of the overdensity threshold. 
Other, true departures from universality may be sought in 
the other components of the universe or in more extreme 
cosmological models. Future extremely well resolved simu¬ 
lations should allow an percent level estimate of the univer¬ 
sality of the mass function. 
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p(EO) 

a 

P 

A 

200p6 

0.6730 ±0.004 

0.1783 ±0.007 

0.4237 ±0.001 

^vir 

0.7369 ±0.004 

0.2089 ±0.007 

0.3894 ±0.001 

200pc 

0.8286 ±0.005 

0.2776 ±0.009 

0.3335 ±0.001 

500pc 

1.0223 ±0.007 

0.3417 ±0.009 

0.2672 ±0.001 

0 

0 

0 

0 

1.2576±0.009 

0.3803 ±0.012 

0.2181 ±0.001 

2000pc 

1.6088±0.015 

0.3824±0.018 

0.1755 ±0.001 


Table Al. Parameters of the best-fitting mass function for EO 
haloes at z = 0. 


APPENDIX A: RESULTS FOR ELLIPSOIDAL 
HALOES 

We identihed haloes using the Ellipsoidal Halo Einder of 
(Despali et al. 2013) and repeated the analysis described in 
the main text. Although the best-ht parameters are system¬ 
atically different, the results are otherwise consistent with 
those found for SO haloes: the virial overdensity yields uni¬ 
versality, and others do not. 

Eigure Al shows the measured mass function for EO 
virial haloes and the residuals with respect to the best ht 
relation (calculated at z = 0), similarly to what was done 
in Eigure 2. Table Al summarises how the best ht z = 0 
parameters depend on overdensity threshold. They behave 
regularly, just as for the SO haloes shown in Eigure 10, and 
are well modelled by equations (Al) (equivalent to the blue 
curves in Eigure 10). 

Eor EO haloes we have obtained analogous relations to 
those of equation (13): 

a = 0.7057+ 0.2125x + 0.3268x^ 

p = 0.2206 +0.0.1937x- 0.04570x^ (Al) 

Ao = 0.3953 -0.1768X. 


APPENDIX B: ISSUES FROM THE INITIAL 
CONDITIONS 

The initial power spectrum P{k) plays an important role 
when translating from halo mass M to the scaled variable 
12 = 6clcr‘^{M)^ since it determines the value of the mass 
variance S{M) = 

In principle, a(M) should be the same for any simula¬ 
tion run with the same cosmological parameters. In practice, 
the number of Eourier modes that can be effectively sampled 
in the initial conditions depends on some computational lim¬ 
its: (i) the box size determines the minimum mode in the 
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Figure Al. Same as Figure 2, but for the haloes in the Ellipsoidal Overdensity catalogue. 


power spectrum that can be sampled by a specific simula¬ 
tion; {ii) each of our simulations started from a different 
random realisation of the displacement field. If not prop¬ 
erly accounted for, these effects may increase the scatter - 
and possibly bias - the halo mass function measured in the 
simulation with respect to theoretical model predictions. 

E.g., using the theoretical linear power spectrum when 
computing a yields a precision which cannot be reduced be¬ 
low about 10%. However, future wide field missions (e.g. 
Euclid Laureijs et al. 2011) require percent level precision. 
As we describe below, to achieve this, we calculate cr(M) 
from the actual realization of the initial power spectrum in 
each simulation box - scaled using linear theory to z = 0 - 
and not from the theoretical linear power spectrum. 


realisations of the displacement field. However, these sim¬ 
ulations also differ in box size, etc. To isolate the effect of 
the random seed, we reran two of the simulations of the sec¬ 
ondary set (“uno” and “wmap7”) using two different random 
seeds. We found few percent-level differences in the high 
mass tail, with one of the two seeds generating greater de¬ 
partures from the universality. See, e.g. the last blue triangle 
in the z = 0 panel of Figure 2 or the gray and brown ones 
in the corresponding panel of Figure 7. This is despite the 
fact that we use the actual realisation of the power spectrum 
when computing a{M). Thus, even when all other parame¬ 
ters are held fixed, some of the scatter in the measured halo 
mass function is introduced by the initial seed. 


B1 Power spectrum of each realisation 

Figure B1 shows some detailed results on the computational 
effects introduced in the initial P{k) regarding (z) the ran¬ 
dom number seed choice and (ii) the box size. The top panel 
shows S — a^(M) as a function of mass M, for all the simu¬ 
lations of the main set; the black line shows the theoretical 
power spectrum from GAME. To compute S(M) for each 
simulation we integrated the initial power spectrum from 
the minimum mode resolvable in each box (see equation (3). 
The bottom panel shows the relative difference between the 
measured S(M) and the one calculated using the theoreti¬ 
cal P(k). For comparison the lower panel shows the initial 
power spectrum for each simulation of our main set and the 
relative differences with respect to the theoretical one. The 
use of the actual power spectrum for each simulation allows 
us to achieve a more precise estimate of ly for the box, and 
hence greater precision on the measured mass function. 


B2 The random seed for the initial displacement 
field 

The six simulations from the main set are independent, 
meaning that they have all been generated from different 


B3 Mass variance definition 


The top part of Figure B2 compares different definitions for 
the mass variance S(M). The relative differences are pre¬ 
sented with respect to the theoretical prediction - i.e. inte¬ 
grating the theoretical power spectrum ( black solid curve). 
The green line in the top panel (orange line in the mid¬ 
dle panel) shows the prediction for S{M) in our smallest 
(largest) simulation Ada (Flora) where we sum over the 
modes in the initial conditions, starting from the k-mode 
corresponding to the box-size. The blue short dashed line 
represents one way of accounting for the limited box-size 
of the simulation: we subtract from the theoretical one the 
mass variance computed using a top-hat filter with scale 


Rbox 


f 3 Mbox 

V 47rp6 



(Bl) 


The red long dashed curve shows the case in which we 
compute S{M) accounting for the fact that the periodic 
boundary conditions of the simulation mean that the Fourier 
modes within the box are constrained to give the background 
density on the scale of the box. This constraint modifies the 
expression for the variance to 
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Figure Bl. Top: Measured mass variance and residuals with 
respect to the theoretical calculation. Bottom: Initial power spec¬ 
trum measured from each simulation and relative difference with 
respect to the theoretical case. 
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Figure B2. Top: Comparison of the variance as a function of 
mass using the different prescriptions described in the text to 
account for the finite box size, for our smallest and largest sim¬ 
ulation boxes. Bottom: Differences between the measured mass 
functions when the prescriptions shown in the top panel are im¬ 
plemented, and when the box size is effectively infinite. 


where 

S.{M) = w[kR(M)]W[kRho.] (B2) 

represents the cross correlation variance between the la- 
grangian scale of the halo and that of the simulation box 
(e.g. Musso & Sheth 2014a,b). For large boxes Sbox s so 
Sx/Shox ~ 1 and the correction to s becomes vanishingly 
small. 

The bottom panels show the z = 0 mass function resid¬ 
uals in Ada and Flora which result from using different pre¬ 
scriptions for accounting for the finite box size (with respect 
to the case in which the box size is infinite). For the larger 


box, Flora, the effect is negligible, but for the smaller box, 
Ada, it is not. If the mass variance does not properly ac¬ 
count for the limited box-size, the rescaling to 12 results in 
non-negligible bias at small masses and large scatter at large 
masses. 

B4 ZA or 2LPT? 

For all our simulations, we used N-GenIC to generate the ini¬ 
tial conditions from a glass distribution; N-GenIC uses first 
order perturbation theory, calculated through the Zel’dovich 
Approximation. Not including second order perturbations in 
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Figure B3. Difference in the measured mass function for different 
IC methods. For four redshifts (in different colours), we show 
the logarithmic residuals between the results obtained using the 
Zel’dovich Approximation (ZA - obtained with N-GenIC) and 
2nd order lagrangian perturbation theory (2LPT - obtained with 
2LPTic) to generate the Initial Conditions. The solid curves are 
calculated from Reed et al. (2013), where the authors present a 
fit to the expected ratio between the two methods, which is in 
good agreement with our measurements. 

the IC has a strong impact on the results obtained at out¬ 
puts “close” to the initial conditions Crocce et al. (2006): 
the differences are larger for simulations with small boxes 
and/or late starting redshifts. We chose z = 99 as a starting 
redshift for all our simulations with a box size larger than 
100 h~^Mpc] for Ada, which has a box of only 62.5 h~^Mpc, 
we generated the ICs at z = 124. To test the impact of ICs 
on our results we ran an identical copy of Bice and of the 
two WMAP7 simulations, using exactly the same parame¬ 
ters and seeds, but generating the IC with 2LPTic instead of 
N-GenIC. The symbols in Figure B3 show the differences in 
the measured mass functions; different coloured solid curves 
show the relation of Reed et al. (2013): 

- 0.12 — 

dnzA/dn2LPT = e , (B3) 

even though it was only calibrated for 1 ^ ^ 5. Note that 

their ly is defined differently than ours: their i^i = yd/ = 
6c(z)/a(M, z). At the four redshifts shown, the differences 
are on average very small - less than 5% at small and in¬ 
termediate masses and less than 10 % in the very high-iy 
tail - which lies within the intrinsic scatter (e.g. that due to 
the random initial seed) seen in the mass function. A small 
difference between the two definitions can be seen at high 
redshift (i.e. at z = 5 in the Figure), but this also lies within 
the intrinsic scatter of the mass function. Since our best 
fits are calibrated using data at z ^ 1.25, they are not af¬ 
fected by second order effects in the initial condition density 
field. For example, the best fit parameters obtained using 
Reed-rescaled points (all z, all cosmologies) are a = 0.7603, 
p = 0.2549 and Aq = 0.330 - similar at per-mil level to the 
ones calculated with the original data. 


APPENDIX C: CLUSTER MASS FUNCTION 

One of the most important application that necessite of 
an accurate and well calibrated mass function is the study 
of the observed cluster counts (Vikhlinin et al. 2009; Rozo 


et al. 2010; Planck Collaboration et al. 2013a; Sartoris et al. 
2015). Degeneracies between fitted parameters mean that 
the best fit parameters to fits restricted to cluster mass ha¬ 
los (M ^ 10^3h~^M q) may differ from those returned from 
fitting a larger range of masses (we refer to these as the 
CMF and HMF, for cluster and halo mass functions, re¬ 
spectively). This Appendix discusses the expected level of 
systematic bias this may induce on cosmological constraints 
derived from the observed number of clusters per square de¬ 
gree. 

The black circles in the top panel of Figure Cl show the 
halo mass function extracted from all cosmologies and red¬ 
shifts from our simulation suite. The main text shows that 
a good fit to the bins that contain at least 30 halos (red tri¬ 
angles) is given by a Sheth & Tormen (1999) mass function 
with the following parameters: Aq = 0.3295, p = 0.2536 and 
a = 0.7689. However it is worth mentioning that accord¬ 
ing to the points density distribution in the considered area, 
this curve may be less accurate in describing the shape of the 
most massive haloes: the cluster mass function. To test the 
difference, and to better describe the objects more massive 
than Mvir ^ 3 x Mq that are typically associated 

with groups and clusters of galaxies, we have performed a 
fit to the orange crosses only. The best fit parameters are 
Aq = 0.8199, a = 0.3141, and p = 0, as reported in Table 3. 
The value p = 0 is due to the fact that it is mainly the small 
masses which determine p (Sheth & Tormen 1999). In the 
first bottom subpanel we show the relative difference - in 
log space - of the two fits for the halo mass function tail. 
For comparison, the other two subpanels show the relative 
residuals of the data points with no error bars to the two 
corresponding fits. 

A relative difference in the number of clusters per square 
degree may appear using the fitting function computed us¬ 
ing only the clusters (hereafter Cluster Mass Function) or all 
haloes over a broader range of masses (Halo Mass Function). 
To quantify this difference in Fig. C2 we present the relative 
difference in cluster counts - Mvir > 3 x - per 

square degree in the Q^m-crs plane between the CMF and the 
HMF fits. To compute the count N we have integrated the 
HMF and CMF fits over Mvir > 3 x IO^^Mq/Zi and comov¬ 
ing volume. The figure shows that in our reference Planck 13 
cosmology the relative difference in the cluster counts be¬ 
tween the two fits is ~ 4%. It increases to ~ 10% for small 
values of both Qrn and as . At high Qrn and as the difference 
between the HMF and CMF-derived counts is smaller. 

In Fig. C3 we present the relative cosmological con¬ 
straints for the cluster counts in the Q^rn-crs plane between 
the HMF and the CMF. The relative difference is well below 
10% in the whole range and the diagonal shape shows the 
degeneration region in the parameter space. 

This suggests that, notwithstanding the small difference 
between the two fits, caution may be necessary when using 
fits to the halo mass function for precision cosmology. For 
cluster counts, the CMF reported in the last row of Table 3 
may be more appropriate than the HMF, and may yield 
better than 4% accuracy. 
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V 


Figure Cl. Halo mass function for all considered cosmologies 
and redshifts extracted from our simulation suite (black circles). 
Red triangles show bins with at least 30 counts while the or¬ 
ange crosses also require that the halos be more massive than 
3 X 10^^ Mq//i. The solid and dashed curves represent the best 
fit to the red triangles and orange crosses - accounting for the 
associated Poisson error bars. The first subpanel shows the rela¬ 
tive difference - in log space - of the two fits, while the other two 
present the relative difference of the data points with respect to 
their corresponding best fit curve. 
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